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' INTRODUCTION 

This  final  report  summarizes  a study  of  hybrid  electro- 
static-electromagnetic techniques  for  the  calculation  of  surface 
currents  caused  by  the  motion  of  photoelectrons  emitted  from 
satellites.  Section  II  presents  the  results  using  the  hybrid 
SQUID  code  for  a series  of  SGEMP  test  cases.  Section  III 
describes  a technique  to  treat  currents  in  small  struts  and 
presents  some  calculations  using  this  technique.  Section  IV 
presents  the  theory  of  nonref lective  boundary  conditions  for 
v electromagnetic  SGEMP  codes.  Included  are  one-  and  two- 

dimensional  calculations  showing  the  effectiveness  of  the 
• boundary  conditions  in  preventing  the  reflection  of  outgoing 

electromagnetic  waves. 
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SECTION  II 


SGEMP  RESPONSE  CALCULATIONS 
USING  THE  SQUID  HYBRID  EM  CODE 

A series  of  SGEMP  response  calculations  was  performed 
using  the  hybrid  SQUID  code  on  the  configuration  shown  in 
Figure  1.  The  diameter  of  the  support  varied  between  zero 
and  10  cm.  This  simple  axisymmetric  configuration  is  similar 
to  that  of  a proposed  SGEMP  experiment.  Our  purpose  was  to 
examine  the  efficiency  of  a hybrid  calculational  technique, 
SQUID,  for  determining  SGEMP  response.  The  SQUID  technique 
treats  heavily  space  charge  limited  behavior  close  to  emitting 
surfaces  with  high  resolution.  It  calculates  dipole  correct 
electromagnetic  response  for  far  field  effects  while  including 
correct  monopole  electric  fields  and  particle  trajectories  near 
emitting  surfaces. 

In  Section  2.1  of  this  report  we  discuss,  in  detail,  the 
SQUID  theory.  The  particular  calculations  performed  are  des- 
cribed in  Section  2.2.  Finally,  we  discuss  the  numerical  re- 
sults and  their  implications  to  future  SGEMP  response  calcula- 
tions in  Section  2.3. 

2.1  THEORETICAL  BASIS  OF  THE  HYBRID  SQUID  CODE 

2.1.1  The  Purpose  of  SQUID 

The  reader  should  be  familiar  with  the  basic  techniques 
of  fully  electromagnetic  (EM)  particle  codes,  which  have  been 
fully  described  elsewhere  (References  1,  2).  This  discussion 


1.  Boris,  J.  P.,  "Relativistic  Plasma  Simulation-Optimization 
of  a Hybrid  Code,"  Proc.  4th  Conf.  on  Numerical  Simulation 
of  Plasmas,  pp.  3-67,  NRL,  Washington,  D.C. , November  1970. 

2.  Katz,  Ira,  James  Harvey  and  Andrew  Wilson,  "Particle  Simu- 
lation Techniques  for  SGEMP,"  Systems,  Science  and  Software 
Report  SSS-R-75-2604 , July  1975. 


Figure  1 - Rxposure  object  and  points  of  interest 
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emphasizes  the  differences  between  SQUID  techniques  and  the 
usual  particle  code  algorithms.  It  is  in  the  neighborhood 
of  object  boundaries  that  rapid  variations  in  field  strengths 
and  charge  densities  occur.  SQUID  is  a technique  for  accurately 
calculating  these  small  wavelength  fluctuations  and  their  effect 
on  the  generation  of  long  wavelength  modes. 

The  SQUID  concept  has  been  developed  to  handle  problems 
in  which  the  spatial  scale  of  electrostatic  fields  close  to 
electron-emitting  surfaces  is  small.  Whenever  the  character- 
istic space-charge  limiting  distance  is  much  less  than  the 
dimensions  of  the  emitting  surface,  an  ordinary  EM  particle 
code  is  incapable  of  treating  the  problem  accurately  unless 
very  small  spatial  zones  are  used,  and  this,  of  course,  leads 
to  very  long  calculations. 

4 

The  initial  ideas  of  the  SQUID  concept  are  outlined  in 
Reference  2.  Briefly  stated,  a one-dimensional  subgrid  is 
defined  in  the  region  along  the  surface  of  a conducting  object 
where  space-charge  limiting  takes  place.  The  particle  dynamics 
in  this  finely  resolved  grid  can  then  be  calculated  accurately. 

The  resulting  charge  density  is  used  to  determine  the  screened 
normal  component  of  the  electric  field  near  the  surface  of  the 
emitting  object. 

Because  the  surface  charge  density  and  normal  electric 
fields  play  an  important  role  in  the  SQUID  method,  the  center- 
ing scheme  adopted  in  the  code  is  different  from  that  of  ordi- 
nary EM  particle  codes.  These  differences  and  the  treatment  of 
boundary  conditions  are  described  in  Section  2.1.2. 

Some  care  must  be  exercised  in  distributing  the  charge 
of  a particle  among  adjacent  cells  in  order  that  the  charge 
dipoles  set  up  along  boundaries  by  moving  particles  are  accu- 
rately calculated.  Currents  on  the  surface  can  be  generated  by 
EM  wave  interactions.  In  Section  2.1.3,  we  describe  an  algorithm 
to  calculate  particle  currents  so  that  transverse  EM  waves  gene- 
rated by  these  individual  particles  will  be  calculated  correct 
to  dipole  order. 
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In  Section  2.1.4,  we  describe  how  physically  meaningful 
longitudinal  electric  fields  can  be  calculated  with  subgrid 
resolution  in  those  cases  where  the  transverse  fields  have 
wavelengths  long  compared  with  cell  dimensions,  but  where  sig- 
nificant charge  density  fluctuations  occur  near  the  surface  on 
a scale  length  small  compared  with  cell  dimensions. 

2.1.2  SQUID  Centering  Scheme 

Conducting  boundary  conditions  on  the  electromagnetic 
fields  in  particle  codes  are  easily  applied;  the  tangential 
electric  field  is  forced  to  be  zero.  In  the  ordinary  centering 
of  particle  code  grids,  conducting  physical  boundaries  pass 
through  points  where  the  tangential  electric  fields  are  defined 
as  in  Figure  2.  Since  the  boundary  conditions  are  applied 
directly,  this  formulation  is  extremely  convenient  to  apply  and 
is  most  commonly  used.  However,  some  quantities  relevant  to 
surface  phenomena,  such  as  surface  charges,  are  not  directly 
accessible  and,  indeed,  must  be  extrapolated  from  interior 
values.  In  dealing  with  the  effects  of  photon  generated  elec- 
trons on  metallic  objects,  the  surface  response  is  the  desired 
information.  Consequently,  we  have  adapted  a modified  center- 
ing scheme  as  shown  in  Figure  3 where  normal  electric  fields 
and  tangential  magnetic  fields  are  defined  directly  on  the  sur- 
face. This  enables  one  to  calculate  surface  charges  and  skin 
currents  directly.  The  boundary  condition  on  the  tangential 
electric  field  is 

E = -E 

II  (below  surface)  ||  (above  surface) 

Since  all  our  field  equations  are  advanced  using  a first-order 
symmetric  difference  form,  this  boundary  condition  is  entirely 
equivalent  to 

e _ i r 

II  (surface)  2 ^E/|  (below  surface)  E||  (above  surface)^ 
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2.1.3  SQUID  Dipole  Correct  Current  Algorithm 
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For  purposes  of  numerical  smoothing,  finite-sized  par- 
ticles are  frequently  used  in  EM  particle  codes.  However,  the 
literal  interpretation  of  a particle  as  a charge  smeared  over 
a finite  region  is  fraught  with  danger  in  problems  where  bound- 
aries play  an  important  role.  In  cases  where  limiting  occurs 
within  a single  cell,  this  interpretation  can  lead  to  incorrect 
results  when  the  dipole  field  (resulting  from  the  motion  of 
charge  away  from  a surface)  is  being  calculated. 

Instead  of  regarding  a particle  as  a finite-sized  object 
which  is  pushed  according  to  appropriate  area  weights  of  the 
electric  fields  in  its  vicinity,  we  prefer  to  consider  the 
particle  as  located  at  a point  in  space  and  to  push  it  by  the 
interpolated  field  at  that  point.  The  difference  in  interpre- 
tation is  merely  semantic  except  near  boundaries.  There,  the 
normal  area-weighted  treatment  must  be  reconsidered. 

While  a particle  is  considered  localized  at  a point  for 
dynamical  purposes,  its  charge  density  must  be  interpolated 
between  adjacent  cells  because  the  charge  density  is  evaluated 
at  specified  points  in  the  grid,  namely,  at  cell  centers  (see 
Figure  3).  Near  boundaries,  special  care  has  to  be  given  to 
this  interpolation  in  order  that  the  current  calculated  gives 
rise  to  the  correct  electromagnetic  field. 

If  a charge  q moves  a perpendicular  distance  x from 
a conducting  surface,  a distant  observer  sees  a dipole  field 
corresponding  to  a dipole  of  strength 

d = 2 qx 

The  modified  centering  scheme  in  Figure  3 corresponds  to  cell 
centered  charge  densities.  First  examine  the  case  of  a 
particle  at  x between  L/2  and  3L/2  from  the  surface  as 
shown  in  Figure  4.  If  this  is  the  only  charge  in  the  system, 
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then  two  equations  are  needed  to  get  the  correct  dipole  strength; 


q=Pl  + P2 


2 qx  = 2 (L/2)  + 2 P2<3L/2) 


Solving  for  and  p2  / find 


’l"1*  (l- 


- L/2 


/x  - L/2\ 
p2  = q \ L / 


Thus,  the  correct  finite  difference  current  j is  just 

!qv  L/2  < x < 3L/2 

0 all  other  x 

For  0 < x < L/2,  there  are  not  enough  unknowns  (p2  = 0)  and  in 
order  to  satisfy  the  distant  observer,  retain  only  the  dipole 


equation 


getting 


2 qx  = 2 P1  (L/2) 


0 < x < L/2 


all  other  x 
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One  can  visualize  these  dipole  correct  currents  as  resulting 
from  finite  size  CIC  particles  of  charge  q centered  at  x 
along  with  image  particles  of  charge  q centered  at  -x.  As 
it  nears  the  surface,  both  the  particle  and  its  image  contribute 
to  the  surface  normal  current  Dgi*  This  enables  one  to  easily 
calculate  currents  for  any  velocity  orientation  which  will  cor- 
respond to  the  correct  dipole  moment  to  a distant  observer. 

The  variation  in  this  dipole  strength  generates  transverse 
electromagnetic  waves. 

2.1.4  Correcting  the  Local  Longitudinal  Field 

In  the  previous  section,  we  dealt  with  the  distant  field 
strength  which  for  a single  emitted  charge  is  dipolar  in  nature. 
However,  the  dominant  field  felt  by  the  test  charge  is  the 
monopole  attraction  from  its  image.  In  the  one-dimensional 
case,  the  field  on  a test  sheet  is  constant  once  the  charge 
leaves  the  surface  and  the  field  on  the  surface  rises  immedi- 
ately to  the  same  constant  as  shown  in  Figure  5.  However, 
the  current  algorithm  described  in  Section  2.1.3  gives  a sur- 
face field  that  rises  linearly  until  the  charge  sheet  is 
half  a cell  out  from  the  surface.  This  is  clearly  an  unphvsi- 
cal  result  when  describing  subgrid  variations  in  charge  density. 
Below,  a technique  is  discussed  for  constructing  the  appropriate 
longitudinal  electric  field  in  the  vicinity  of  the  surface 
assuming  that  the  charge  density  varies  much  more  slowly 
parallel  to  the  surface  than  in  the  normal  direction. 

The  charge  density,  a , on  a surface  zone  can  be  expressed 
as  the  sum  of  time  integrals  of  the  skin  current  k and  perpen- 
dicular particle  current  j , 


°--f 


dt  (\7,.k  + V(.  j) 
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Figure  5 - The  surface  electric  field  for 
a charge  sheet  a distance  x 
above  the  conducting  surface. 
The  solid  line  is  the  physical 
field,  the  dashed  line  is  that 
generated  using  the  dipole  con 
sistent  current  described  in 
Section  2.1. 


If  the  dipole  correct  current  discussed  in  Section  2.1  is  used, 
the  surface  charge  reflects  both  the  linear  turn-on  of  the  sur- 
face electric  field  and  the  partial  tangential  current  cancel- 
ling due  to  the  image  current. 

<7,.  . = -/  (V  • k,.  . +V,  • j,.  . ) dt 

dipole  J ||  dipole  _L  dipole 


It  is  simple  to  calculate  in  the  code  the  actual  particle  cur- 
rent leaving  the  surface.  The  image  charge  contribution  to 
the  surface  current  may  be  calculated  in  the  following  manner. 
If  a particle  of  charge  q is  a distance  0 < x < L/2  from 
the  surface,  its  contribution  to  the  parallel  current  (real)  is 


= qv 


n 


However,  for  the  dipole  current,  the  image  contribution  is  inclu- 
ded in  the  current  and 


.dipole 

3il 


particle  image 


Therefore,  in  order  to  find  the  true  (monopole)  charge  density 
on  a surface,  include  the  complete  image  particle  current 


(image)  * q^l| 

in  the  skin  current  k.  Assuming  that  the  contributions 
to  the  surface  charge  from  particles  more  than  one  cell 
away  are  the  same  in  both  monopole  and  dipole  approximations, 
then 
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( 


°monopole  °dipole  + 


Jmonopole  •'dipole^ 


+ V,  • (k  , - k . ) 1 dt 

li  monopole  dipole  J 


Knowing  amonopole'  true  surf^ce  electric  field  can  be 

found.  Spatial  variation  of  the  monopole  electric  field  normal 
to  the  surface  can  easily  be  approximated  if 


I 9x  1 1 1 * 9x^  I 

by  a simple  integral 


E,(x)  = E±(0)+  / p(5)d4 


This  equation  is  the  basis  of  the  SQUID  subgrid  screening  ap- 
proximation. 

In  practice,  a pseudo-one-dimensional  charge  density  is 
calculated  for  cells  adjacent  to  emitting  surfaces  as  shown  in 
Figure  6.  These  charge  densities  are  then  used  to  give  a 
screened  normal  component  cf  the  electric  field  which  acceler- 
ates the  particles. 


El(x.)  = e|  f + 
1 _L  surface 


f ' 


(x^ ) dx^ 


t+1/2  t— 1/2  A1(  , v ^ t 

v,,  =vM  + At(e/m)  E . 


t+1/2  t-1/2  . f/  , . -t,  , 

v,  = v,  + At(e/m)  c (x, ) 

i 1 1 
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The  Electromagnetic 
Grid 


SQUID  Grid 


Spacecraft 

Surface 


The  grid  for  the  high  fluence  SGEMP  problem.  Each 
cell  adjacent  to  a boundary,  n,  contains  an  entire 
quasi-one-dimensional  grid.  Associated  with  each 
SQUID  grid  is  a one-dimensional  charge  density 
cQ(xj_)  used  for  calculating  the  normal  component 
or  the  electric  field  within  the  cell  adjacent 
to  the  soacecraft  surface. 


As  a result,  the  electrostatic  field  in  the  limiting  sheath 
region  is  resolved  in  great  detail  in  the  direction  normal 
to  the  surface,  the  direction  in  which  it  varies  most  rapidly. 


2.2  HYBRID  EM  CALCULATIONS 

Calculations  of  the  response  of  an  isolated  cylinder  to 
photogenerated  electron  emission  have  been  performed  in  accor- 
dance with  Task  1 (Section  4.1)  of  the  Work  Statement  using 
the  SQUID  code. 

For  the  sequence  of  calculations,  a regular  mesh  with 
equal  spacings  of  10  cm  in  the  r and  z directions  was  em- 
ployed for  electromagnetic  purposes.  The  five  normal  emission 
zones  were  subdivided,  in  accordance  with  the  SQUID  formulation, 
into  0.5-cm  regions.  The  appearance  of  the  gridded  problem  is 
shown  in  Figure  7. 

Electron  emission  was  accomplished  through  the  use  of 
seven  energy  bins.  The  angular  distribution  f(0,4>)  was 
cos9  in  the  normal  direction  and  uniform  in  the  polar  angle 
<p.  The  bin  structure  is  as  follows: 


Bin 

Fractional 
f luence 

Energy  (keV) 

Velocity  (in  units  of 
10  keV  velocity) 

1 

0.20028 

2.79 

0.50758 

2 

0.24054 

7.43 

0.86174 

3 

0.21810 

12.44 

1.1139 

4 

0.16546 

17.34 

1.3157 

5 

0.96678 

22.21 

1.4896 

6 

0.7181 

29.79 

1.7227 

7 

0.00713 

45.94 

2.1398 

The  bin  selection  was  strictly 

random  for  each  particle 

in  most 

runs,  but  in 

the  final  series. 

conditional  selection 

was  employed  to  guarantee  that  an  equal  number  from  each  bin 


i 

j 


i 
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4 


was  chosen  each  cycle.  The  emission  «was  uniformly  spaced  in 

r along  the  surface.  The  pulse  shape  was  given  by  a triangle 

-9 

with  rise  time  of  2.5  x 10  sec  with  a peak  current  of 
2 

7.88  amp /cm  . 

The  major  parameter  variations  among  the  completed  runs 
have  been  the  number  of  particles  emitted  per  cycle,  the  method 
of  bin  selection  and  the  time  step.  Benchmark  results  were 
obtained  using 

At  = 5 x 10  ^ sec 

emission  = 300  particles/cycle 

This  calculation  had,  on  occasions,  over  6000  particles  within 
the  grid,  and  the  200  cycles  were  required  to  go  out  to  10  ns. 
Total  Central  Processor  Unit  (CPU)  time  for  the  run,  and  a 
myriad  of  diagnostics,  was  1622  sec.  Subsequent  runs  were 
primarily  to  determine  the  accuracy  of  reducing  the  numbers 
of  particles  emitted  and  increasing  the  time  step.  Since  the 
time  step  is  not  Courant  limited  by  the  electromagnetic  calcu- 
lation as  long  as  it  is  below  2.2  x 10  ^ sec,  a 10  l0  sec  time 
step  was  chosen  as  an  optimum  from  both  cost  and  pulse  resolution 
viewpoints.  By  reducing  the  particle  emission  to  70  per  cycle 
and  using  the  10  sec  time  step,  CPU  time  dropped  almost  an 
order  of  magnitude  to  225  sec,  with  an  accuracy  of  a few  per- 
cent in  quantities  characterized  by  long  wavelengths.  This  was 
achieved  by  using  nonindependent  energy  bin  selection,  that  is, 
by  constraining  the  energy  bin  choice  so  that  an  equal  number 
from  each  bin  is  emitted  per  cycle.  Surface  currents  in  the 
cylinder  midplane  are  shown  in  Figure  8.  The  upper  curve  is 
from  the  benchmark  run,  the  lower  from  the  few  particle,  long 
time  step  run.  The  closeness  of  these  results  implies  that,  at 
least  for  the  simple  geometry  considered,  accurate  SGEMP  surface 
currents  can  be  obtained  with  very  short  computations. 
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Figure  8 - Midplane  surface  currents  as  a function  of  time.  The 
upper  curve  is  the  benchmark  calculation;  the  lower  u 
conditional  bin  selection  and  required  only  3 minutes 
CPU  time. 


The  production  runs  had 


At  = 10  ^ sec 

emission  = 140  particles/cycle 

Typical  runs  used  5 minutes  of  CPU  time  and  had  a maximum  of 
around  1500  particles  in  the  mesh  at  any  one  time.  The  maximum 
deviation  from  the  desired  net  emission  current  was  less  than 
6 percent,  while  for  70  percent  of  the  cycles,  the  emission  cur- 
rent was  within  3 percent  of  the  desired  triangular  ramp. 

The  5 cm  radius  support  rod  was  inserted  by  setting  the 

electric  field  component  in  the  z-direction  to  zero  for  the  first 

10  zones  near  the  axis.  For  the  1 cm  radius  support,  a set  of 

radial  electric  fields  and  6 component  magnetic  fields  1 cm 

off  axis  was  constructed.  Maxwell's  equations  were  solved 

for  these  positions  using  first-order  differences.  However, 

the  [ (1/r) (3/3r) rB , ] contribution  to  E near  the  axis 

p z 

was  no  longer  a centered  first-order  difference  and,  as  a 

result,  was  accurate  only  to  order  AR  (as  opposed  to  the  rest 

of  the  grid  where  the  centered  differences  were  accurate  to 
2 

order  AR  ) . 

Skin  currents  at  points  C and  E (see  Figure  1)  are 
shown  in  Figures  9 and  10.  These  currents  are  in  amperes 
per  meter  and  are  obtained  by  multiplying  the  surface  magnetic 
fields  (in  gauss)  by  79.53.  The  most  prominent  feature  of  the 
midplane  currents  (Figure  9)  is  the  second  peak  occurring 
approximately  5 nanoseconds  after  the  first  peak.  One  possible 
mode  that  would  give  use  to  this  second  peak  is  that  of  a wave 
generated  at  point  B,  going  out  to  the  wall  radially  and  being 
reflected  back  in.  The  path  length  of  such  a wave  differs  by 
about  1.5  meters  from  the  direct  B to  C distance. 
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In  Figure  10,  the  inclusion  of  a conducting  support 
allows  charge  to  escape,  preventing  the  formation  of  a "sloshing" 
mode.  As  the  diameter  of  the  support  increases,  its  inductance 
decreases,  so  that  the  current  transmitted  reaches  higher  values. 

Complete  time-dependent  information  for  points  A through 
G has  been  plotted  and  is  available  upon  request. 


2.3 


DISCUSSION  OF  RESULTS 


In  this  subsection,  two  aspects  of  the  computational  tech- 
niques used  in  performing  these  calculations  are  examined:  (1) 

Is  subgrid  screening  important  in  determining  particle  dynamics 
and  SGEMP  response  variables  (such  as  skin  currents)?  (2)  Is 
the  hybrid  approach  an  efficient  technique,  or  are  ordinary  EM 
codes  competitive  when  it  comes  to  computer  costs? 

The  importance  of  subgrid  screening  can  be  shown  in 
several  ways.  Analytically,  it  is  an  easy  procedure  to  show 
that  if,  for  a monoenergetic  electron  beam  in  an  EM  code, 

Ax  is  the  EM  grid  spacing,  v is  the 
emission  velocity  and 


Ax  >>  v/co  where 
P 

to  the  plasma  frequency  of  the  beam, 

P 

the  electric  field  within  the  first  zone  will  be  described  in- 


correctly since  the  limiting  distance  is  less  than  a zone  size. 
Since  the  finite  difference  representation  of  Maxwell's  equa- 
tions assumes  linear  variation  of  electric  fields  within  a zone, 
the  finite  difference  representation  is  clearly  inappropriate. 
For  the  case  in  point,  more  than  half  the  current  is  composed 

9 

of  electrons  below  12  keV  in  energy  (v  % 7 x 10  cm/sec) . 

2 

An  emitted  current  of  1 A/cm  (1/7  of  peak  value)  gives 
a plasma  frequency 

o ir,9  -1 
co  a 2 x 10  sec 
P 
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and  a limiting  distance 


v/u)  ~ 3.5  cm 

P 

, Therefore,  the  majority  of  the  electrons  should  limit 

well  within  the  10  cm  zone  size  for  almost  the  entire  duration 
of  the  pulse.  Consequently,  one  would  expect  significant  field 
variations  within  zones  bordering  the  emitting  surface. 

Figure  11  is  a plot  of  the  normal  component  of  the  electric 
field  on  a particle  in  the  near  axis  zone  as  a function 
of  distance  from  the  emitting  surface.  The  departure  from 
linearity  is  obvious.  Figure  12  shows  the  positions  of  the 
particles  at  3 nanoseconds  into  the  calculation.  The  high 
degree  of  limiting  within  10  centimeters  is  apparent. 

How  this  nonlinear  electric  field  behavior  reflects  it- 
self in  measurable  quantities  was  examined  by  running  the  zero 
radius  support  calculation  without  the  subgrid  resolution.  As 
is  expected,  the  induced  currents  are  much  greater,  and  are 
almost  a factor  of  3 more  than  calculated  with  SQUID  grid. 

The  dipole  field  is  incorrect  near  the  surface  since  it  exerts 
too  little  force  on  low-energy  particles,  and  as  a result,  their 
contribution  to  the  dipole  current  is  unphysically  large.  We 
conclude  that  if  large  electromagnetic  zones  are  used  in  this 
calculation,  subgrid  electric  field  corrections  must  be  made. 

Even  though  the  SQUID  technique  enables  one  to  employ 
large  electromagnetic  zones,  the  particle  pushing  and  field 
calculations  are  more  complicated  than  conventional  EM  codes. 
The  relevant  question  then  concerns  the  cost  effectivenss  of 
the  hybrid  technique.  In  this  area,  the  results  seem  quite 
encouraging . 

Based  on  the  spatial  resolution  required  to  perform  an 
accurate  treatment  of  the  heavily  space-charge-limited  region 
near  the  emitting  surface,  the  SQUID  zone  chosen  was  0.5  cm. 
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The  perpendicular  component  of  the  electric  field  as  a function  of 
distance  from  the  surface  3 nanoseconds  into  the  calculation.  An 
magnetic  grid  zone  is  10  cm  long.  The  field  strength  varies  in  no 
fashion  an  order  of  magnitude  within  the  first  electromagnetic  cod 
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Figure  12  - The  electron  sheath  3 nanoseconds  into  the 
problem.  A large  grid  cell  is  10  cm,  the 
SQUID  grid  spacing  is  0.5  cm. 


- . ■ — ZZL*.. . . ■ -- 


The  EM  grid  spacing  was  10  cm.  In  order  to  achieve  the  same 
resolution  by  a purely  "brute-force"  means  it  would  have  been 
necessary  to  reduce  the  EM  grid  spacing  and  hence  the  time 
step  by  a factor  of  20.  In  fact,  the  savings  obtained  by  the 
SQUID  technqiue  were  less  than  this  on  account  of  the  additional 
complexity  in  pushing  particles  in  the  SQUID  region.  However, 
this  increased  the  cycle  time  by  less  than  a factor  of  2,  thus 
resulting  in  a net  savings  of  an  order  of  magnitude. 

It  is  worth  considering  the  economy  of  the  SQUID  calcu- 
lations performed  here.  Runs  were  made  with  time  steps  of 
5 x 10  ^ sec  and  10  sec  and  the  number  of  particles  emitted 
per  cycle  ranging  from  60  to  300.  The  few  particle  (70)  runs 
using  conditional  energy  bin  selection  and  the  longer  10  sec 
time  step  gave  results  within  10  percent  of  those  from  the  long- 
est run  on  long-wavelength  variables  (e.g.  , induced  skin  cur- 
rents) while  requiring  less  than  5 minutes  of  UNIVAC  1108  time 
to  run  including  long  diagnostic  outputs  and  plots . This  trans- 
lates to  running  costs  of  the  order  of  $15.  The  significance 
of  the  short  running  time  is  that  three-dimensional  calculations 
performed  in  an  analogous  fashion  would  be  extremely  practical 
in  contrast  with  the  computer  time  and  storage  requirements  of 
three-dimensional  calculations  using  ordinary  EM  particle  codes. 

We  conclude  that  for  these  highly  space-charge-limited 
SGEMP  calculations  the  SQUID  hybrid  technique  has  been  extremely 
effective.  The  results  strongly  indicate  that  the  dipole-SQUlD 
theory,  as  discussed  in  Section  2.1  of  this  report,  contains 
the  physics  relevant  to  SGEMP.  While  a very  fine  zoned  electro- 
magnetic particle  code  also  contains  the  relevant  physics,  it 
automatically  follows  unimportant  short-wavelength  transverse 
fields.  The  time  step  requirements  of  SQUID  are  of  the  order  6 
of  the  lowest  energy  electrons  times  that  of  an  ordinary  EM  code 
of  comparable  spatial  resolution.  For  the  examples  presented, 
all  of  which  have  0.5  cm  resolution,  the  conservative  10  ^ sec 
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time  step  is  more  than  a factor  of  6 larger  than  the  largest 
stable  time  step  for  an  ordinary  electromagnetic  code.  Thus, 
SQUID  enables  one  to  save  a large  percentage  of  computer  time 
with  no  loss  in  accuracy. 


31 


SECTION  III 


A TECHNIQUE  FOR  TREATING  SMALL  CURRENT 
CARRYING  STRUTS  IN  ELECTROMAGNETIC  PARTICLE  CODES 


The  prediction  of  the  electromagnetic  response  of  a 

satellite  to  an  external  impulse,  caused  for  example  by  a 
— 8 

short  (10  sec)  pulse  of  X rays,  presents  a formidable 
computational  problem.  The  state-of-the-art  for  the  self- 
consistent  determination  of  trajectories  of  electrons  and 
the  electromagnetic  field  generated  by  them  is  well  developed 
for  many  two-dimensional  problems.  Where  there  are  highly 
disparate  geometrical  length  scales,  however,  as  in  a hollow 
right  circular  cylinder  of  radius  R with  an  axially  centered 
strut  of  radius  a <<  R , either  the  computational  mesh  must 
be  extremely  fine  or  other  methods  for  treating  small  scale 
lengths  must  be  developed.  In  this  section,  methods  are 
formulated  for  treating  thin  struts  (a  in  the  above  examples) 
in  both  two-  and  three-dimensions.  A primary  result  of  the 
formulation  is  the  retention  of  a spatial  difference  mesh 
that  is  appropriate  for  the  larger  scale  length,  together 
with  the  relatively  greater  economy  in  computation  associated 
with  the  longer  timestep  that  is  permitted  (by  the  Courant 
stability  condition)  in  the  staggered-leapf rog  finite  dif- 
ference scheme  for  solving  Maxwell's  equations. 

For  the  purposes  which  follow,  the  Maxwell  equations 
governing  the  evolution  in  time  of  the  electric  and  magnetic 

_ — V — V 

fields  E and  B,  respectively,  are  written  in  cgs  units 

1 3B 
' c 3t 

ill  1 SE 

c c 3 1 

where  j is  the  current  density  and  c the  velocity  of  light. 


7 * E = 
V x B = 
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3.1 


FIELDS  NEAR  A STRUT 


First  consider  the  fields  in  the  vicinity  of  a per- 
fectly conducting  strut  of  radius  a . In  a cylindrical 
coordinate  system  whose  z-axis  concides  with  the  strut  axis, 


3E 
3 z 


3E 


Near  the  strut. 


(PfZ,t) 


E (p,z,t)  = 
2iTaBg  (a,  z , t) 


3 1 


aB0 (a , z , t) 

P 

aEp  (a , z , t) 

P 

4ttI  (z  ,t) 
c 


(1) 

(2) 

(3) 

(4) 


and 


2iTaE  (a,z,t)  = 4irQ(z,t)  (5) 

P 

The  current  I carried  by  the  strut  and  its  charge  Q per 
unit  length  are  related  by  the  continuity  equation 


9Q  . *SI 
3 1 3 z 


0 


(6) 


Here,  we  have  assumed  that  the  strut  neither  emits  nor  absorbs 
charge . 

Integrating  Equation  (1)  from  a to  p gives 
P P 

V»  ■ h f Vc',d0'  + Ik  / (7> 

a a 
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since  Ez(a)  = 0*  From  Equations  (2)  through  (5) 


Ez(p)  = 2in  (p/a)  |S  + 'L-  An  (p/a)  |i 


(8) 


The  area  average  of  E z over  a cell  is  given  by 


E = c'1  ifi  + L £1 

z,av  ^ 3z  L 3t 


(9) 


where 


C-1  = 2 <£nfp/a)> 


(10) 


and 

l = y c-1  (id 

are  the  capacitance  and  inductance,  respectively,  per  unit 
length  of  the  strut.  Here,  <f>  denotes  an  area  average  of 
f over  a cross  section  of  the  cell  normal  to  the  strut.  If 
the  cross  section  of  area  A is  approximated  by  a cylindrical 
cell  of  radius  R = centered  on  the  strut  axis,  then 

R 

C_1  = JP  J P^n(p/a)dp  (12) 

a 

= 2 An(R/a)  + a2/R2  - 1 (13) 

In  obtaining  Equation  (9) , it  has  been  assumed  that  the 
dominant  local  fields  are  quasi-static  and,  moreover,  are  axi- 
symmetric  in  the  coordinate  system  centered  on  the  strut.  These 
assumptions  exclude  very  short  wavelength  modes  corresponding, 
for  example,  to  variations  depending  on  the  azimuthal  angle. 

Such  short  wavelength  excitations  are  probably  not  very  impor- 
tant in  any  practically  occurring  SGEMP  context. 
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3.2  3-D  DIFFERENCE  EQUATIONS 

For  three-dimensional  (3-D)  calculations,  we  adopt 
Cartesian  geometry  and  a computational  grid  whose  cells  are 
identical  and  whose  three  mutually  perpendicular  edges  have 
arbitrary  lengths.  The  centering  scheme  is  depicted  in  Figure 
13.  Letting  (i,j,k)  be  the  coordinates  of  a cell  center  in 
units  of  the  lengths  of  the  basic  cell  edges,  we  center  the 
field  and  current  components  according  to  the  following  scheme: 


, E : 
X 

(i+!s , j ,k+J*) 

■V 

(i » j+J5 ,k+^) 

■V 

(i,  j ,k) 

B : 
X 

(i,  j+%,k) 

B : 

y 

(i+35 , j ,k) 

B : 
z 

(i+% , j+% ,k+%) 

Using  a staggered-leapfrog  scheme  in  which  E and  B 
are  advanced  in  a fully  space-and-time  centered  manner,  the 
difference  forms  of  Maxwell's  equations  become 


(i  , j ,k)  - Bfc  ^ (i  , j ,k)  = -cAt  (VpXE) t (i  , j ,k) 


(14) 


Et+1(i,j,k)  - Et(i,j,k)  = cAt  (VpXB) t+i*  (i  , j ,k) 

4tt  tt+*s  , . . . > 

- — 3 (i,3 ,k) 


(15) 


The  subscript  D denotes  the  difference  gradient, 


(VDxA)x(i, j ,k)  = 


Az  (i,  j+*5,k)  -Az  (i,  j-*s,k) 

_ 


Ay(i,j  ,k+*s)-Ay  (i,  j ,k-h) 
Az 


(16) 
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(VDxA)y (i, j ,k) 


A (i,  j,k+*s)-A  (i,  j ,k-h) 

A A 

Az 


A (i+h, j,k)-A  (i-%, j,k) 

Z Z 

Ax 


(17) 


(vdxa)z  (i, j,k) 


Ax  (i' j+Js»k)_Ax  (i,  j-^/k) 
— 


Ay  ( i+ht j f k ) “Ay (i-h, j , k) 
Ax 


(18) 


and  i,j,k  are  integral  or  half-odd  inteqral  according  to  the 
centering  scheme  of  Figure  13. 

The  foregoing  equations  are  intended  to  apply  for 
systems  where  there  is  no  great  disparity  in  length  scales. 

For  systems  with  highly  dissimilar  length  scales,  such  as  the 
cylinder  with  struts,  the  foregoing  equations  must  be  modified 
to  avoid  the  fine-zoning  associated  with  a small  strut  radius. 
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3.3  LINEAR  STRUT  IN  A 3-D  GEOMETRY 


Consider  a linear  strut  whose  axis,  parallel  to  the 
z-direction,  pierces  the  x-y  plane  at  the  position  (I,J)  where 
Ez  is  centered.  The  grid  is  sufficient  to  resolve  spatial 
variations  of  the  electromagnetic  field  except  within  the  zones 
through  which  the  strut  passes.  There,  the  local  fields  vary 
according  to  Equations  (3)  and  (4).  Thus,  the  difference  Equa- 
tions (14)  and  (15)  apply  at  all  space  points  and  for  all  field 
components  with  the  following  exceptions: 

a)  the  equations  for  Bt+^  (I , J±35,k) -Bt  2 (I , J±^,k) 

X X 

and  By+h  (IiJj,  J,k) -B^  (I±^,J,k) 

and 


4-  4-  "I  +* 

b)  the  equations  for  E^  (I , J ,k) -E^ (I , J ,k)  . 


The  required  modifications  of  the  equations  for  advancing  the 
magnetic  field  components  B^  and  B^  are  obtained  by  setting 
E (I,J,k)=E  in  Equation  (14).  Thus,  for  example. 


B^+'2  (I , J +h , k ) -B^  ’-(I,J+b,k)  = 

X X 


r E*(l,J+l,k)-E*  (I , J ,k) 

' ' civ 


e£  (I , J+»s,k+Js)  -Ey  (i , J+>s,k-Js) 
Az 

(19) 


Now  it  remains  only  to  alter  the  equations  for  advancing  E (I,J,k) 
in  time.  The  required  modifications  are  based  upon  Equations  (6) 
and  (9)  together  with  Ecuation  (20) , 
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1 Ezlav(I-J-l°-El;,avII"J'k) 

c At 


- P-  J^U.J.k)  - p-  1 


t+a 


(k) 


+ (VxB)^(I,J,k) 


(20) 


where  E is,  again,  the  area  average  of  E over  a cross 

Z / cl  V Z 

section  of  the  zone  normal  to  the  strut. 


Tt+a  Tt+1  , 

I = ctl  + (l-a)i 


is  the  current  carried  by  the  strut, and  the  parameter  a indi- 
cates the  degree  of  implicitness.  For  a = 1(0),  the  current 
is  treated  in  a fully  implicit  (explicit)  manner;  for  a = h, 
the  current  is  time  centered  and  the  treatment  of  the  strut 
current  is  half-implicit.  The  two  additional  Equations,  (6) 
and  (9) , are  sufficient  to  describe  the  system  which  now  con- 
tains the  two  additional  variables  I(z)  and  Q(z). 

For  advancing  I and  Q,  Equations  (6)  and  (9)  are 
written  in  the  fully  space  and  time  centered  form 


T It+1(k)-It(k)  . 

L At  C 


Qt+h (k+h) -Qt+h (k-h) 

+ \ 

Et+1  (k) +Efc  (k) 

Az 

2 

z , av  z , av 

(21) 


(k+^ ) -Qt~^  (k+i;  1 


At 


Ifc  (k+l)-It (k) 
Az 


(22) 


Here,  the  dependence  on  all  discretized  spatial  variables  except 
those  indicating  position  along  the  strut  has  been  suppressed. 

The  stability  of  the  computational  scheme  corresponding 
to  Equations  (20)  through  (22)  with  a = 1 and  a = 4 has  been 
demonstrated  in  the  two-dimensional  calculations  described  in 
Section  6.  In  the  context  of  stability,  it  is  worth  noting 
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that  Equations  (20)  through  (22)  depart  from  the  explicit 
staggered-leapf rog  scheme  that  is  commonly  employed  in  the  solu- 
tion of  Maxwell's  equations.  If  the  conventional  scheme  were 
retained  throughout,  then  I and  Q would  be  centered  at 
t + and  t,  respectively.  Equation  (20)  would  have  a = h 
and  Equations  (21)  and  (22)  could  be  replaced  by  the  explicit 
forms 


. It+Js(k)-It_!5(k) 

L At 

Qt+1  (k+*s)  -Qfc  (k+h) 

At 

For  sufficiently  small  values  of  L,  i.e.,  for  rods  of  radius 
comparable  with  the  cell  size,  we  conjecture  that,  in  contrast 
to  Equations  (21)  and  (22) , the  stability  of  these  equations 
require  a time  increment  more  restrictive  than  that  imposed  by 
the  Courant  conditions  over  the  portion  of  the  mesh  through 
which  struts  do  not  pass;  viz, 


= -C-1 


Qt  (k+^-Q1  (k-*s) 


Az 


+ E (k) 
z ,av 


Az 


(23) 
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3.4  TWO-DIMENSIONAL  SIMULATIONS  OF  STRUT  GEOMETRY 


Since  off-axis  struts  in  an  otherwise  axisyimnetric 
geometry  pose  an  inherently  three-dimensional  problem,  the 
development  in  the  previous  section  was  three-dimensional  from 
the  outset.  On  the  other  hand,  in  a cylinder  with  identical 
struts  parallel  to  the  z-axis  located  at  the  same  radial  dis- 
tance from  the  cylinder  axis,  and  equally  spaced  in  angle, 
the  azimuthally  averaged  fields  for  the  mode  of  highest  angular 
symmetry  satisfy  precisely  the  two-dimensional  field  equations 
of  axisymmetric  cylindrical  r,z  geometry.  The  mode  of  highest 
angular  symmetry  is  the  only  one  excited,  for  example,  by  a 
radially  uniform  pulse  of  X rays  incident  normally  on  one  end 
of  the  cylinder. 

Particular  interest  is  attached  to  the  first  two  of  the 
relevant  Maxwell  equations  given  below  in  the  vicinity  of  the 
radial  location  of  the  struts: 


1 

c 


i 3B.  9E 

9E 

1 0 r 

z 

c 9 1 9 z 

9r 

9E„  . 

_ 5.  + iH  j 

9 1 c Jz 

1 9 

r 9r 

<rV 


1 

c 


+ 


^0 
9 z 


(25) 

(26) 
(27) 


Here,  the  indicated  fields  are  azimuthally  averaged  and  r is 
measured  from  the  axis  of  the  cylinder  (not  from  a strut  axis) . 

The  difference  equations  in  the  strut-bearing  zones,  which 
have  radial  index,  I,  are 
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Bt+'^  (l±h , i+h)  (i±jj,1+is) 

At 
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1 ^+1(I,j)-E^(I,j)  + p.  + p 3 (I/j) 


At 


c 2irrIAr  c Jz 


1 (rI+  j)-(rI-  |^)Bt+Js(I-Jsf  j) 
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(29) 


The  quantity  I ( j ) is  the  total  current  in  all  of  the  struts 
at  the  z location  with  index  j . 

To  complete  the  description  requires  a determination  of 
the  current  I;  the  required  relations  will  be  given  in  the 
form  of  Equations  (6)  and  (9) . Equation  (6)  is  simply  the  equa- 
tion of  continuity  along  the  strut.  Equation  (9)  must  be  derived 
from  a heretofore  unused  Maxwell  equation  which  describes  depar- 
tures from  azimuthal  symmetry.  The  presence  of  current  and 
charge  on  the  struts  lead,  in  fact,  to  a radial  magnetic  field, 
B^,  which  in  the  highest  symmetry  case  considered  here,  vanishes 
on  averaging  over  azimuth.  The  defining  equation 
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when  integrated  twice  over  0 along  the  curve  r = constant  on 
which  the  strut  axes  lie,  yields 
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where  N is  the  number  of  struts  and  0 the  angle  subtended 

o 

at  r = 0 by  a strut  radius. 

For  struts  separated  by  a distance  b <<  r,  it  is  reason- 
able to  approximate  and  Eg  in  the  manner  described  in 

Section  3.1,  in  which  case  Br  and  Eg  are  proportional  to  the 
current  and  charge  per  unit  length,  respectively,  on  the  struts. 
This  is  a basically  quasi-static  approximation  and  requires, 
as  well,  that  the  wavelengths  of  the  fields  being  computed  are 
also  large  compared  to  the  strut  separation.  In  practice,  this 
is  not  always  the  case,  but  the  considered  approximations  are 
about  as  well  as  we  can  do  in  the  two-dimensional  simulation 
of  an  inherently  three-dimensional  problem.  Results  based  on 
these  two-dimensional  approximations  are  presented  in  Section 
3.5. 


3.5 


TEST  CALCULATIONS 


In  order  to  test  the  difference  formulations  derived  in 
the  previous  sections,  a series  of  axially  symmetric  2-D 
electromagnetic  calculations  was  performed.  They  fall  into 
two  general  classes.  The  first  consisted  of  tests  of  "loaded" 
struts,  that  is,  small  diameter  struts  connecting  large  objects 
so  that  the  dominant  features  of  the  system  are  the  strut  induc- 
tance and  the  objects'  capacitance.  The  second  tests  involved 
the  behavior  of  an  isolated  radiating  rod.  Here,  the  system 
is  inherently  electromagnetic  in  nature  with  the  observable 
features  being  the  frequency  of  the  radiation  and  the  radiation 
"resistance"  damping  oscillations  on  the  rod. 

For  the  "loaded"  strut  tests,  a series  of  runs  was  made 
of  a conducting  cylinder  supported  by  off-axis  cylindrically 
symmetric  struts  in  a conducting  tank.  The  configuration  is 
shown  in  Figure  14. 

The  dimensions  of  the  system  are  the  following: 

Inner  cylinder:  radius  83.5  cm 

length  100  cm 

Outer  cylinder:  radius  2.5  m 

length  6 m 

Support  struts:  length  2.6  m 

distance  off-axis  41.75  cm 

The  calculations  were  performed  using  FRED/EM-PRECHARGE , a 
2-D  fully  electromagnetic  code  which  has  the  capability  of 
initializing  the  fields  to  satisfy  Laplace's  equation.  The 
calculational  mesh  was  15  x 30  with  AR  = 16.7  cm  and 
AZ  = 20  cm,  and  the  timestep  chosen  was  3.85  x 10  second, 
which  is  0.9  of  the  Courant  timestep.  Initially,  the  object 
was  charged  to  10  kV  with  respect  to  the  tank  using  the 
PRECHARGE  routine,  an  SOR  Laplace  solver  that  generates  elec- 
tric fields  which  are  curl  free  to  machine  accuracy.  The 
dominant  capacitance  of  this  configuration  is  just 
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6 00  cm 


Figure  14  - The  configuration  used  for 

testing  the  response  of  loaded 
struts . 
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1.31  x 10  Coulomb 

2 

10  volts 


1.3  x 10 


-10 


farads 


This  corresponds  to  a capacitance  of  117  cm  (esu) . 

The  charge,  q,  on  the  interior  cylinder  was  calculated  by 
summing  the  normal  (and  only)  component  of  its  surface  electric 
field.  This  is  done  automatically  in  the  PRECHARGE  routine. 


For  a first  test,  the  strut  inductance  per  unit  length 
-21  2 2 

was  set  at  5 x 10  sec  /cm  . This  corresponds  to  a quasi- 
static period  of 


= 2tt/lC  = 2ttX 


260  x 5 x 


10  21) 


117  sec  = 7.75  x 10 
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sec 


or  201  code  cycles.  We  would  expect  a priori  the  actual  period 

to  be  somewhat  longer  since  the  light  transit  time  across  the 

— 8 

tank  is  2 x 10  sec.  The  particular  value  of  L was  chosen 
to  give  a period  long  enough  to  be  clearly  resolved  by  the 
code,  but  not  so  long  that  the  calculations  became  ponderous. 
Using  Equations  (11)  and  (13)  , we  find  that  this  corresponds 

to  a case  of  a single  2.38  cm  diameter  strut.  The  actual  code 

*"8  — 8 

period  was  9.34  x 10  sec  or  1.6  x 10  sec  longer  than  the 
quasi-static  prediction.  This  is  quite  plausible  considering 
that  the  period  required  for  information  to  pass  back  and  forth 

_ O 

through  the  tank  is  4 x 10  second.  The  current  midway  along 
the  strut,  as  a function  of  time,  is  shown  in  Figure  15.  The 
initial  conditions  excited  several  modes,  accounting  for  the 
nonsinusoidal  behavior.  That  the  electric  field  obtained  from 
Equation  (20)  is  indeed  appropriate  comes  from  examining  the 
magnetic  field  a half  zone  outside  the  strut  at  a radius  of 

_ O 

50.1  cm.  For  example,  at  6.11  x 10  sec,  the  current  in  the 
strut  80  cm  from  the  tank  wall  was  69.6  amperes.  Cylindr ically 
averaged,  this  would  generate  an  average  static  magnetic  field 
of  0.2777  gauss  at  50.1  cm.  The  code  shows  a field  of  0.2786 
gauss,  extremely  close  to  the  static  value.  Since  only  the 
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Figure  15  - Strut  current  vs. time  for  an  induction  per  unit  length  of 
5 x 10~  sec  2/cm  . 


free  space  Faraday's  equation  is  used  to  generate  the  magnetic 
field,  it  shows  clearly  the  correctness  of  the  electric  field 
chosen  at  the  strut.  This  kind  of  correlation  accuracy  is 
maintained  throughout  the  calculation.  To  minimize  the  dif- 
ference, the  magnetic  field  used  as  an  example  was  chosen  when 
the  displacement  current  was  very  small. 

As  another  example,  the  oscillations  of  a case  with  an 
inductance  one  sixteenth  of  the  first  case  were  calculated. 

This  corresponds  to  a symmetrical  support  of  eight  struts, 
each  1.6  cm  in  diameter.  Since  the  predicted  quasi-static 
period  is  smaller  than  a transit  time,  the  calculated  dynamic 

electromagnetic  period  is  much  longer  than  the  predicted  quasi- 
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static  period,  7.3  x 10  sec  vs.  1.9  x 10  sec.  Also,  the 
free  space  inductance  (not  just  the  magnetic  field  energy  in 
the  strut  zone  itself)  is  comparable  to  the  strut  inductance 
which  would  lengthen  the  period  in  the  electromagnetic  calcu- 
lation. Strut  current  vs.  time  from  this  calculation  is  shewn 
in  Figure  16 . 

A simpler  formulation  of  the  "loaded"  strut  equations, 

neglecting  the  capacitance  of  the  strut,  was  tested.  If  the 

results  were  similar  to  those  of  the  previous  calculations, 

neglecting  the  capacitance  would  be  a useful  approximation 

since  it  requires  less  computer  storage.  However,  there  were 

serious  difficulties  using  the  inductance  alone.  In  the  case 
-21  2 2 

of  L = 5 x 10  sec  /cm  , the  period  was  lengthened  by 
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0.7  x 10  sec  (Figure  17).  More  importantly,  the  current  on 
the  strut  was  not  uniform  as  a function  of  position  along  the 
strut.  This  resulted  in  charge  buildup  on  the  strut  and  un- 
physical electric  fields.  Figures  18  and  19  shew  electric  field 
plots  both  without  and  with  the  strut  capacitance.  Notice  the 
large  fields  near  the  strut  in  the  no  capacitance  case  (Figure 
18)  and  the  very  small  fields  near  the  strut  in  the  capacitance 

case.  The  plots  were  selected  in  both  cases  to  maximize 
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fields  near  the  struts.  The  absence  of  the  capacitance  pre- 
vents current  from  flowing  with  a uniform  velocity  down  the  strut. 

The  radiative  properties  of  an  isolated  strut  in  a tank 
were  also  tested  using  the  FRED/EM  code.  The  first  configur- 
ation is  shown  in  Figure  20.  The  discretized  strut  included 
six  charges  and  five  currents.  The  zone  si^es  were  as  before. 
Initially,  the  charge  was  distributed  sinusoidally  on  the  strut. 

Current  vs.  time  at  the  strut  midpoint  is  shown  in  Figure  21. 

-9 

The  period  is  approximately  8 x 10  sec,  corresponding  to  a 
rod  length  of  120  cm.  Figure  22  shews  the  same  configuration 
with  the  radial  zone  size  increased  to  30  cm.  This  increases 
the  distance  to  the  radial  wall  resulting  in  a "clear"  time 
of  3 x 10  sec  as  compared  to  1.7  x 10  sec  in  the  case  of 

Figure  21.  Figure  23  shows  a configuration  with  the  larger  £R 

and  a rod  increased  in  length  by  two  cells  (40  cm) . The  period 

-9 

is  now  10.8  x 10  sec,  corresponding  to  a strut  length  of 
162  cm.  These  times  have  an  uncertainty  of  order  1/2  of  time- 
step  (which  was  3 x 10  ^ sec).  This  translates  to  4.5  cm 
uncertainties  in  the  strut  length.  Thus,  within  the  uncertain- 
ties, the  first  two  cases  radiated  like  a 120  cm  (6  zone)  strut, 
and  the  second  like  a 160  cm  (8  zone)  strut.  The  nonsinusoidal 
behavior  is  due  to  the  initial  conditions  exciting  more  than 
the  lowest  mode. 

One  feature  worth  noting  in  these  isolated  strut  calcu- 
lations is  the  decay  in  amplitude  of  the  strut  current  up  until 
the  reflected  signal  returns.  This  is  completely  expected  due 

to  the  power  lost  from  radiation.  For  the  2.38  cm  diameter, 

_ 8 

120  cm  strut,  theory  predicts  a damping  time  of  1.3  x 10 
second.  Estimates  of  the  actual  numerical  damping  rate  based 
on  Figure  22  agree  well  with  the  theory.  However,  more  defini- 
tive calculations  are  required  to  accurately  determine  the 
numerical  damping  period. 
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Figure  22  - Current  on  a 120  cm  isolated  strut.  Tank  ra 


Time  (sec) 

on  an  isolated  160  cm  strut.  Tank  radius  is 


3 . 6 CONCLUSIONS 
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We  have  presented  both  theory  and  calculations  support- 
ing a technique  for  simultaneously  solving  free  space  wave 
equations  along  with  quasi-static  circuit  elements.  The  cal- 
culational  results  show  that  the  technique  is  both  stable  and 
accurate.  The  simplicity  of  the  final  equations  (e.g.,  (21) 
and  (22) ) makes  them  very  easy  to  incorporate  in  existing 
Maxwell's  equation  codes.  Since  the  new  equations  appear  to 
have  little  effect  on  timestep  size,  current  flowing  through 
small  cross  section  objects  can  be  treated  without  resorting 
to  fine  zoning  and  its  concurrent  timestep  limitations.  Thus, 
in  case  2,  presented  in  Section  3.5  where  the  strut  diameter 
was  1.6  cm,  we  were  able  to  use  a timestep  of  3.845  x 10  ^ sec, 
more  than  a factor  of  7 larger  than  would  be  allowed  with 
fine  zoning.  Equally  important,  the  formulation  presented 
allows  the  two-dimensional  calculation  of  the  m = 0 response 
of  three-dimensional  systems,  such  as  off-axis  symmetrical 
struts.  Included  properly  are  effects  of  field  penetration 
between  regions  exterior  and  interior  to  the  radial  location 
of  the  struts. 

We  feel,  therefore,  that  these  techniques  make  practical 
the  computation  of  SGEMP  response  on  many  new  classes  of  systems, 
in  particular,  those  systems  with  current-carrying  cables  or 
struts  that  have  cross  sectional  dimensions  small  compared  to 
electromagnetic  wavelengths. 
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SECTION  IV 


NONREFLECTIVE  FREE  SPACE  BOUNDARY  CONDITIONS 
FOR  SGEMP  CODES 


The  grid  in  SGEMP  codes  is  only  of  finite  extent.  For 
most  SGEMP  calculations  that  have  been  carried  out,  the  light 
transit  time  across  the  grid  is  short;  that  is,  comparable 
to  the  driving  pulsewidth.  If  waves  are  back-scattered  or 

reflected  from  the  edge  of  the  grid,  they  can  interact  with 
fields  and  currents  in  the  emitting  region.  For  a satellite 
in  free  space,  the  "clear  time"  for  reflections  should  be  long 
compared  to  the  characteristic  times  associated  with  the 
driving  pulse  so  that  back-scattered  waves  do  not  enter  into 
the  calculated  response.  So,  in  order  to  calculate  SGEMP 
response  properly,  it  is  necessary  to  devise  boundary  conditions 
which  do  not  reflect  electromagnetic  waves. 

One  technique  which  has  been  used  is  to  impose  a multi- 
pole outgoing  wave  boundary  condition.  That  is,  the  fields 
on  the  outermost  edge  are  varied  as  if  the  fields  in  their 
neighboring  cells  were  outgoing  spherical  waves.  In  general 
practice,  only  the  lowest  multipoles  are  included.  While, 
conceptually,  this  technique  is  appropriate,  in  actual  use 
it  has  some  flaws. 

One  difficulty  is  that  in  order  to  use  only  the  low 
order  multipole  terms,  the  boundaries  should  be  a large  distance 
(compared  to  object  dimensions)  from  the  radiation  source.  This 
requires  either  a large  number  of  zones  or  large  zones  nearer 
the  boundary.  A second  difficulty  is  that  the  waves  generated 
numerically  do  not  all  have  phase  and  group  velocity,  c,  which 
is  independent  of  wave  number.  This  is  especially  true  for 
wavelengths  comparable  to  a few  times  the  grid  spacing.  As  a 
result,  especially  if  the  grid  spacing  is  made  large  near  the 
outer  boundary,  there  will  be  phase  mismatching  on  the  boundary. 
A third  difficulty  concerns  the  cumulative  effects  of  errors. 
Since  the  procedure  is  to  generate  fields  at  the  boundary  which 
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exactly  cancel  reflected  wave  amplitudes,  any  phase  error 
results  in  both  the  reflected  wave  not  being  cancelled  as 
well  as  a new  incoming  wave  being  generated.  As  a result, 
those  modes  whose  reflections  are  not  eliminated  well  can  be 
pumped  by  the  boundary  conditions  leading  possibly  to  unstable 
numerics . 

Our  approach  to  the  nonref lecting  boundary  problem  has 
been  to  try  amplitude  damping  as  opposed  to  phase  matching. 

The  reasoning  behind  this  is  rather  transparent;  even  if  one 
does  not  damp  a wave  completely,  at  least  one  can  avoid  pump- 
ing the  mode . 

If  the  region  of  space  is  current  free,  then  the  two 
curl  equations  of  Maxwell's  equations  can  be  written 

3E 

yjr  = + c V x B (30) 


3B 

= - c V X E (31) 


Adding  a damping  coefficient,  v,  to  each  equation,  we  get 


3E 

= c V x B -vE  (32) 


3B 

= - c 7 x E -vB  (33) 

O L 


If  the  spatial  dependence  of  v is  assumed  to  be  extremely 
small,  then,  if  EQ(x,t)  and  BQ(x,t)  are  solutions  to 
Equations  (30)  and  (31),  the  solutions  to  Equations  (32)  and 
(33)  can  be  written  as 
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wr 


-vt 


E(x,t)  = EQ(x,t)  e 


B(x,t)  = Bo(x,t)  e“Vt 


As  a result,  waves  in  such  a damped  region  would  decay  in 
time . 

Finite  difference  equations  corresponding  to  Equations 
(32)  and  (33)  are 


Et+1  - Et  = At(c  V x B - vE)t+1/2 


= cAt  V x Bt+1/2  - vAt 


tet+1  + Et 


,t+l 


(l  - ^■)Et  + cAt  V x Bt+1/2 
vAt\ 


(1  + 


2 / 


B 


t+1/2 


(l  _ li|tjBt_1//2  - cAt  V x Et 

pfj 


(34) 


(35) 


Time  centering  is  preserved  and  a degree  of  implicitness 
is  introduced  which  makes  the  equations  stable.  In  practice, 
v is  a function  of  position  and  is  only  nonzero  very  close  to 
the  outside  of  the  mesh. 

In  the  following  section  are  analytical  results  describing 
the  damping  and  reflection  of  waves  obeying  Equations  (32)  through 


j 

i 
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(34)  for  both  continuum  and  discrete  cases.  We  present  the 
results  of  a series  of  one-dimensional  numerical  calculations 
performed  to  determine  the  optimal  spatial  variation  of  the 
damping  factor,  v, in  Section  4.2.  The  last  section  shows  the 
effectiveness  of  this  technique  in  preventing  reflected  waves 
in  a two-dimensional  calculation. 


4 . 1 ANALYTICAL  RESULTS 

It  is  instructive  to  analyze  the  steady  oscillations 
corresponding  to  Equations  (32)  and  (33)  and  to  their  discrete 
analogues.  Equations  (34)  and  (35),  for  a one-dimensional  infinite 
medium  with  v > 0 in  the  half  space  x > xq  and  vanishing 
in  the  space  x < xq. 

Assuming  a steady  wave  propagating  in  the  positive 
x-direction  for  x > xq,  we  shall  determine  the  relative  ampli- 
tudes of  left  and  right  going  waves  in  x < xq. 

For  the  continuous  case,  the  substitution 


E (x , t ) = s(x,t)cosh  t - B(x,t)sinh  t 
B (x , t ) = 3(x,t)cosh  t - e(x,t)sinh  t 


v (x ' ) dx  ' 


reduces  Equations  (32)  and  (33)  to 


equations  which  are  formally  equivalent  to  the  original  Maxwell 
equations  without  the  damping  term. 


Equations  (39)  admit  a steady  wave  propagating  to  the 


right 
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i (oit-kx) 


3 


e e 
o 

i (o>t-kx) 

e e 


e 


with 


co  = ck 


(40) 


(41) 


Thus , 

E = B = ee-T  (42) 

also  propagates  to  the  right.  The  E and  B fields  are  damped 

in  x > x ; in  x < x , the  fields  have  their  free  space  values 
o o 

and  there  is  no  reflected  wave. 

In  the  discrete  case,  the  waves  may  also  be  damped  in 
x > xq;  but, in  general,  there  will  exist  a reflected  wave  in 
the  half  space  x < xq.  Now  assume  that  v = constant  for 
x > x^,  i.e.,  for  the  discrete  spatial  coordinate  n >_  N, 
and  vanishes  otherwise. 

Assuming  a solution  of  the  form 


E = E e 
o 


B = B e 
o 


i ( ju)A  t-nk  ' Ax) 


i [ ( j + 1/2)  o)At  - (n+1/2)  k ' Ax] 


(43) 


Equations  (34)  and  (35)  reduce  to 


/ 2 i sin 
\ “ 


p)At 

2 , ooAt 

+ v cos  — 


)B°- 


2i  c sin 


k'Ax 


Ax 


(44) 


(- 


sin 


At 


u)At 


+ v cos 


03 At  \ 

2 


E = 


2i  c sin 
Ax 


k ' Ax 


(45) 
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where 


E = ±B  ; in  particular,  for  E = B 
o o o o 


k’Ax 

Ax 

toAt 

_ 1 

vAx 

o)At 

sin  2 

cAt 

sin  2 

-L 

2c 

V—vJ  O 

2 

kAx 

vAx 

wAt 

= 

sin 

— - 1 

2c 

cos 

2 

(46) 

where 


sin 


kAx 

2 


Ax 

cAt 


sin 


wAt 

2 


(47) 


is  the  dispersion  relation  in  the  absorption-free  half  space. 

The  imaginary  part  of  k1  determines  the  attenuation  in 

Ax 

the  absorbing  medium.  Writing  k'^—  = ur  + ^1'  anc^  separating 
real  and  imaginary  parts  in  Equation  (46),  there  results 


cosh  Uj. 

sin 

UR  = 

kAx 
sin  — 

sinh  Uj 

cos 

UR  = 

vAx 

- 2c  COS 

coAt 

2 

By  inspection , 

we  f 

ind 

G 

M 

A 

O 

for 

kAx 

2 

t=  IfN 

vl 

and 


< 1 
Ax 

Thus,  all  wavelengths  resolvable  on  the  mesh  are  damped;  wave- 

c A t 

lengths  near  A = 2Ax  are  only  weakly  damped  for  near 

C t ^ a 

unity  and  are  undamped  for  ~~  = 1. 

To  obtain  the  reflection  coefficient,  consider  solu- 
tions having  the  form  of  Equation  (43)  with  Bq  = for 

n > N,  and  the  form 
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1 


I 


E = ei3“At  (ae~iknAx  + beiknAx) 

n i(j  + l/2)a)At  , -ik(n+l/2)Ax  , ik(n+l/2)Ax> 

B = e J (ae  - be  ) 


for  n <_  N-l . 

Equations  (43)  determine  E^  and  Bn+1/2’  Tke  comPlex 
coefficients  a and  b are  determined  using  Equation  (34)  for 
n = N and  Equation  (35)  for  n = N-l,  remembering  that 
vN  = v and  vN_^  = 0.  Following  a simple  but  tedious  calcula- 
tion, we  find  the  reflection  coefficient 


r = 


b 

a 


-i (k'-k) 

1 - e 

-i (k'+k) 

1 + e 


Ax 

2 


Ax 

2 


where  the  complex  wave  number  k'  is  related  to  the  given  real 
quantities  k and  v by  Equation  (46) . 

In  general,  although  Equation  (46)  must  be  solved  numeri- 
cally before  r can  be  determined,  two  limiting  cases  are 
amenable  to  simple  analysis.  When 


k 


1 


and 


then 


vAx 

2c 


<<  1 


cAt  „ 
Ax  - 


x 


r 


vAx 
* 2c 


c A t 

For  = 1 and  kAx  = tt,  r = 0.  Nevertheless,  since 

k'  = 0,  a wave  encountering  a boundary  at  finite  x in  the  damp- 
ing medium  will  be  reflected  and  return  with  undamped  amplitude 
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1 


to  the  absorption  free  region.  In  practice,  however,  < 1, 

and  there  will  be  some  damping  of  the  reflected  wave,  the 
amount  depending  on  the  thickness  of  the  absorbing  medium. 
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4.2 


1-D  PARAMETER  STUDY 


The  spatial  variation  of  the  damping  coefficient,  v, 
in  Equations  (32)  and  (33)  makes  the  problem  of  analytically  pre 
dieting  wave  behavior  difficult  since  the  space  is  no  longer 
homogeneous.  As  a preliminary  indication  of  the  method's 
validity,  a one-dimensional  transverse  electromagnetic  wave 
code  was  written.  This  code  has  100  cells  in  the  x-direction 
with  edge  centered  y-components  of  electric  field  and  cell 
centered  z-components  of  the  magnetic  field.  It  numerically 
follows  the  propagation  of  plane  polarized  TEM  modes  in  the 
x-direction.  The  wave  source  was  a sinusoidally  varying 
y-component  current  located  at  the  midpoint  of  the  grid.  The 
damping  coefficient,  v,  was  zero  everywhere  except  within  a 
specified  number  of  cells  from  either  edge  of  the  mesh.  Cal- 
culations were  run  to  determine  the  functional  form  of  the 
damping  coefficient's  spatial  dependence  which  best  minimizes 
the  amplitude  of  the  reflected  wave  after  the  wave  source  was 
turned  off  (50  cycles) . 

After  several  runs,  several  characteristics  of  the 
technique  were  apparent.  The  first  was  that  increasing  the 
number  of  damping  zones  from  5 to  10  on  each  side  of  the  mesh 
made  small  differences  (factors  of  2)  in  the  reflected  wave 
amplitude.  Second,  the  damping  should  increase  slowly  from 
zero  to  some  number  such  that 

vAt  > 1.0 

Finally,  the  damping  effectiveness  was  proportional  to  wave- 
length . 

Let  us  examine  the  results  of  a set  of  runs  made  with 

— j = 3.8.  The  damping  increases  quadratically  from 

' max 

zero  in  each  boundary  region 


v (x)  At 
2 


3.8*  (5-x)  /25,  x <5 
0 , 5 ^ x 4 95 


, u , o $ 

v 3 . 8 (x- 


(x-95)  /25 , x >95 


Table  1 presents  the  maximum  electric  field  amplitude  in 
the  mesh  after  200  timesteps  (cAt  = 0.90)  for  a range  of  pump 
frequencies.  It  is  apparent  that  longer  wavelengths  damp  most 
effectively,  as  discussed  in  Section  4.1,  but  that  even  for  high 
frequency,  short  wavelengths,  over  93  percent  of  the  wave 
amplitude  is  damped  after  200  cycles.  The  damping  of  the 
longest  wavelength  is  almost  99  percent  effective.  Figure  24 
shows  a plot  of  the  percent  remaining  after  200  cycles  vs. the 
frequency  of  the  pump  mode.  The  relation  is  linear  to  a 
remarkable  degree.  This  bodes  extremely  well  for  SGEMP  res- 
ponse calculations,  since  the  wavelengths  associated  with  the 
strong  electromagnetic  waves  generated  by  the  ringing  of  surface 
currents  are  usually  large  compared  to  grid  spacings. 

One  point  to  note  is  that  the  damping  coefficient  gets 
considerably  larger  than  the  timestep,  that  is, 


vAt  > 1 


”rv:;r 


Table  1 

THE  EFFECTIVENESS  OF  THE  DAMPING 
ON  WAVES  FOR  VARIOUS  FREQUENCIES 


1/wAt 

|B  (t=50)  | 

|B(t=200)  1 
1 'max 

% remaining 

10 

0.550 

0.037 

6.7 

15 

0.556 

0.024 

4.3 

20 

0.521 

0.018 

3.5 

30 

0.511 

0.012 

2.3 

60 

0.504 

0 . 006 

1.2 
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4.3  DAMPED  WAVE  BOUNDARIES  IN  TWO-DIMENSIONS 

To  test  the  usefulness  of  the  damping  theory  in  two 
dimensions,  a sample  calculation  was  performed.  Previously, 
the  behavior  of  a precharged  rod  inside  a conducting  cylinder 
had  been  calculated.  The  system  parameters  were 


cylinder 

length  = 

600  cm 

radius  = 

450  cm 

rod 

length  = 

120  cm 

radius  = 

1.19 

mesh 

Az  = 20 

cm 

AR  = 30 

cm 

The  configuration  is  shown  in  Figure  25. 

The  current  at  the  middle  of  the  wire  as  a function  of 
time  is  shown  in  Figure  26.  In  this  case,  the  precharge  was 
chosen  to  be  purely  sinusoidal  as  opposed  to  those  in  Section 
3.5,  and  as  a result,  very  little  harmonic  distortion  is  evident 


Damping  was  added  to  the  outer  time  zones  using  the 
quadratic  prescription  as  in  the  one-dimensional  example.  In 
this  case,  we  calculated  a damping  coefficient  based  on  the 
z-coordinate  and  another  based  on  the  r-coordinate . The  damp- 
ing coefficient  used  was  the  maximum  of  the  two.  Figure  27  show 
how  the  damping  effectively  reduced  the  reflections  from  the 
walls.  High  frequency  distortion  is  apparent  when  the  effect 
of  reflections  was  first  seen  in  Figure  27.  This  results  from 
the  inefficient  damping  of  high  frequency  modes  as  discussed 
in  the  previous  sections.  The  overall  stability  of  the  system 
should  be  noted;  there  are  no  growing  oscillations  at  late  time 
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Figure  26  - The  current  at  the  wire  center  vs.  time  for 

the  case  with  no  damping. 


Strut  current 


Figure  27  - The  current  at  the  wire  center  vs. time  with  damping 

in  the  outermost  5 zones.  The  damping  parameters 
were  chosen  from  the  results  of  the  1-D  study 
described  in  Section  4.2. 
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